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We systematically generalize the exotic ®He-B phase, which not only exhibits unconventional 
symmetry but is also isotropic and topologically non-trivial, to arbitrary partial-wave channels 
with multi-component fermions. The concrete example with four-component fermions is illustrated 
including the isotropic /, p and d-wave pairings in the spin septet, triplet, and quintet channels, 
respectively. The odd partial-wave channel pairings are topologically non-trivial, while pairings in 
even partial-wave channels are topologically trivial. The topological index reaches the largest value 
of in the p-wave channel {N is half of the fermion component number). The surface spectra 
exhibit multiple linear and even high order Dirac cones. Applications to multi-orbital condensed 
matter systems and multi-component ultra-cold large spin fermion systems are discussed. 
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Superconductivity and paired superfluidity of neu¬ 
tral fermions possessing unconventional symmetries are 
among the central topics of condensed matter physics. If 
Cooper pairs formed by spin-i fermions carry non-zero 
spin, their orbital symmetries are in the odt^artial-wave 
channels. The p-wave paired superfluidity [3 includes 
the ^He-A phase exhibiting point nodes Q, the fully 
gapped B phase Qjand the recently reported polar state 
with linear nodes [1, . The p-wave superconductiviW 

has also been extensively investigated in SrRu 2 04 [ 3 - 
0, and heavy fermion systems including UGe 2 , URhGe, 
UCoGe M- The p-wave superfluid ^He and supercon¬ 
ductors exhibit rich topological structures of vortices and 
spin textures under rotations or in external magnetic 
fields, respectively [ll|, [l^- In addition, experimental 
signatures of the possible nodal /-wave superconductiv- 
ity have also been reported in UPta [l^ [l^ . 

Among these unconventional pairing phases, the ^He- 
B phase is distinct: in spite of its non-s-wave pairing 
symmetry and spin structure, the overall pairing struc¬ 
ture remains isotropic and fully gapped. Its pairing ex¬ 
hibits the relative spin-orbit ^mmetry breaking from 
S'Ol(3) ® SOs{5) to SOj{3) [if where L, S, and J rep¬ 
resent the orbital, spin, and total angular momentum, 
respectively. The relative spin-orbit symmetry-breaking 
has also been studied in the context of Pomeranchuk in¬ 
stability termed as unconventional magnetism leading to 
dynamic generation of spin-orbit coupling [Il[i3- 

Furthermore, the ^He-B phase possesses non-trivial 
topological properties [IMl- Topological states of mat¬ 
ter have become a major research focus since the dis¬ 
covery of the integer quantum Hall effect [20l-l^. Re¬ 
cently, the study of topological band structures has ex¬ 
tended from time-reversal (TR) breaking systems to TR 
invariant systems (25 - 1^ . from two to three dimensions 
flM and from insulators to superconductors [13- 

Il9l. bsMSlf . The ^He-B phase is a 3D TR invariant topo¬ 
logical Gooper pairing state. Its bulk Bogoliubov spectra 
are analogous to the 3D gapped Dirac fermions belonging 


to the Dill class characterized by an integer-valued in¬ 
dex [T^. The non-trivial bulk topology gives rise to the 
gapless surface Dirac spectra of the mid-gap Andreev- 
Majorana modes [s^l • Evidence of these low energy states 
has been reported in recent experiments [s^. 

Because the electron Gooper pair can only be either 
spin singlet or triplet, the p-wave ^He-B phase looks the 
only choice of the unconventional 3D isotropic pairing 
state. In this article, we will show that actually there 
are much richer possibilities of this exotic class of pairing 
in all the partial-wave channels of L > 1. We consider 
multi-component fermions in both orbital-active solid 
state systems and ultra-cold atomic systems with large 
spin alkali and alkaline-earth fermions, both of which 
have recently attracted a great deal of attention 
For simplicity, below we introduce an effective spin s to 
describe the multi-component fermion systems with the 
component number expressed as 2N = 2s -I- I > 4. Gom- 
pared with the 2-component case, their Gooper pair spin 
structures are greatly enriched [H, H^. For example, 
the 4-component spin-| systems can support the /-wave 
septet, p-wave triplet, and d-wave quintet pairings, all of 
which are fully gapped and rotationally invariant. Never¬ 
theless, only the odd partial-wave channel ones, i.e., thep 
and /-wave pairings are topologically non-trivial. Their 
topological properties are analyzed both from calculating 
the bulk indices and surface Dirac cones of the Andreev- 
Majorana modes. For the p-wave case, the topological 
indices from all the helicity channels add up leading to a 
large value of N"^. Gorrespondingly the surface spectra 
exhibit the coexistence of 2D Dirac cones of all the orders 
from 1 to 2N — I. 

We begin with an /-wave spin septet Gooper pairing 
Hamiltonian in a 3D isotropic system of spin-| fermions 

H = ^e^J^{k)Co.{k) - ^ Pl^MPm,u{k')A^) 

fr2l2 

in which -p and p is the chemical potential. 

a = ±|, is the spin index, g is the pairing interaction 
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strength, and Vq is the system volume. The pairing oper¬ 
ator is defined as = c]^{k)Y 3 rnik)[S^'"I{\apc^p{—k) 

where k = k/k, Y 3 m{k)’s with —3 < m < 3 are the 3rd 
order spherical harmonic functions, and with —3 < 
V <i are the rank-3 spherical tensors based on the spin 
operator S in the spin |-representation, where v is the 
eigenvalue of Sz- For later convenience, Y^m(k) are nor¬ 
malized according to \ Y 3 m(k)\'^ = 1. i? is the charge 
conjugation matrix defined as Rap = (—)‘^^'^5a-p sat¬ 
isfying RS^R~^ = —S such that Rapc^p transforms in 
the same way under rotation as Ca does. The expres¬ 
sions for spherical harmonic functions and spin tensors 
are presented in Appendix 

After the mean-field decomposition, Eq. [T] becomes 

^ ( 2 ) 

k 

in which k is summed over half of momentum space; 
'k(fc) = (eg is the Nambu spinor; the order pa¬ 

rameter Am,u is defined through the self-consistent equa¬ 
tion as 

A™,. = ^ ^{G\c^{-k)Yi^{k)R^S^''’^cs{k)\G) (3) 

k 

with (G'|...|G) meaning the ground state average. The 
matrix kernel H{k) in Eq. [2] is expressed as 

H{k) = e{k)T 3 0 / 4 x 4 + A(fc)r+ -|- A(-fc)T_, (4) 

where and t± = i(ri ± 1 x 2 ) are the Pauli matrices 
acting in the Nambu space. A(A;) is defined in the matrix 
form in spin space as 

A{k)=J2{S^''R)d*’''{k), ( 5 ) 

where d*'''{k) = Am,vYzm{k) and is dubbed as the d- 
tensor in analogy to the d-vector in ^He. The usual d- 
vector is represented in its three Cartesian components, 
while here, the d-tensor is a rank-3 complex spherical 
tensor. 

We consider the isotropic pairing with total angular 
momentum J = 0, which is a generalization of the p- 
wave ^He-B phase. Similarly, it is fully gapped, and 
thus conceivably energetically favorable within the mean- 
field theory. Its d'^(fc) can be parametrized as d'^(fc) = 
CfAf(-^)^Y 3 ^{k), where c/ is an overall normalization 
factor given in Appendix [BJ Ag is the complex gap mag¬ 
nitude, or, equivalently, 

Aik) = Afi^fKf{k)R ( 6 ) 

in which Kf = CfU{k)S^°W(k); U{k) rotates the z- 
axis to k as defined in the following gauge U{k) = 




FIG. 1: Pictorial representations of the pairing matrices 
over the Fermi surfaces of (a) the isotropic /-wave septet 
pairing and (fo) the isotropic p-wave triplet pairing with 
spin-1 fermions. Intuitively, the /-wave matrix kernels 
U{k)S^^U\k) and the p-wave ones U{k)S^^U^{k) for each 
wavevector k are depicted in their orbital counterpart har¬ 
monic functions in (a) and (fe), respectively. 

g-^<Pks^g-^ekSy which 9k and (j)k are polar and az¬ 
imuthal angles of k, respectively. The explicit form 
of A(fc) and the corresponding spontaneous symmetry 
breaking pattern are presented in Appendices [B] and ICl 
respectively. 

With the help of the helicity operator h{k) = k ■ S, 
Kf{k) can be further expressed in an explicitly rotational 
invariant form as 

KfCk) = -h^Ck) + ^hCk), (7) 

which is diagonalized as W{k)Kf{k)U{k) = (—fS'f + 
^S'z). For a helicity eigenstate with the eigenvalue 
A, the corresponding eigenvalue of Kf(k) reads 
^ respectively. 

The Bogoliubov quasi-particle spectra are E\{k) = 
+ |A/P(A) 6^2 satisfying Ex{k) = E_x{k) due 
to the parity symmetry. 

Next we study the pairing topological structure. The 
pairing Hamiltonian Eq. Uin the Bogoliubov-de Gennes 
(B-deG) formalism possesses the particle-hole symmetry 
GpE[{k)G~^ = —E[*{—k) with Gp = n 0 R. Further¬ 
more, the isotropic pairing state described by Eq. [S] 
is TR invariant satisfying CTH(k)G^^ = H*{—k) with 
Gt = h® R, and thus it belongs to the DIB class. The 
associated topological index is integer-valued which will 
be calculated following the method in Ref. [3l| . H(k) is 
transformed with only two off-diagonal blocks as e{k)Ti -|- 
Af{-^)^K{k)T2. The singular-value-decomposition to 

its up-right block yields U{k)L{k)A{k)W{k), in which 
L{k) and A(fc) are two diagonal matrices only depen¬ 
dent on the magnitude of k defined as Lxx{k) = E\{k) 
and A\x{k) = respectively. The angles satisfy 

ta,n9x{k) = and for simplicity A^: is set as 

positive. The A:^-dependence of the pairing amplitude is 
regularized: Beyond a cutoff k^, Af vanishes. 
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FIG. 2: The gapless surface spectra for the isotropic /-wave 
septet pairing in (a) and for the isotropic p-wave triplet pair¬ 
ing in (b) with spin-| fermions. 

The topological index is calculated through the SU(4) 
matrix = U(k)A(k)U^k) as 

^- = ^1 ( 8 ) 

which is integer-valued characterizing the homotopic 
class of the mapping i.e., 7r3(5'17(4)) = Z. Neverthe¬ 
less, Nyj is only well-defined up to a sign: After changing 
Af —>■ —Af, flips the sign. As shown in Appendix 
m at /r > 0 iVu, is evaluated as 

N^= X! ^ sgn(^A)- (9) 

A=±|,±l 

Its dependence on sgn(^A) is because 9x{k) varies from 
0 —>' 7 rat^A >0 but from tt — 7 > ^ 0 as fc 

varies from 0 to k f to +oo. A similar form of Eq. [^was 
obtained in Ref. pflj in which the Fermi surface Chern 
number plays the role of A in Eq. [9l For two helicity 
pairs of A = ±| and A = ±^, their contributions are 
with opposite signs, and thus = 2. 

The non-trivial bulk topology gives rise to gapless sur¬ 
face Dirac cones. Because of the pairing isotropy, with¬ 
out loss of generality, an open planar boundary is chosen 
a,t z = 0 with /i(z) = e/> 0 atz <0 and —oo at 
z > 0. The mean-field Hamiltonian becomes H(k\\,z) in 

which k\\ = {kx, ky) remains conserved while the transla¬ 
tion symmetry along the z-axis is broken. The symmetry 
on the boundary is C„oo including the uni-axial rotation 
around the z-axis and the reflection with respect to any 
vertical plane. C^ao is also the little group symmetry at 
k\\ = 0, then the four zero Andreev-Majorana modes at 
k\\ = 0 are Sz eigenstates denoted as |0a,/). The associ¬ 
ated creation operators 7 ^ are solved as 



-f e *(2 + I)c_„(fc|| = 0, z)]uq,(z), (10) 

where (p is the phase of A; Uf^a{z) is the zero mode wave- 
function exponentially decaying along the z-axis, and its 
expression is presented in Appendix [El 


The surface zero modes |0 q,,/) at fey = 0 possess an 
important property that the gapped bulk modes do not 
have: They are chiral eigen-modes satisfying Cch\0a,f) = 
| 0 q-j) with Va = 0 for a = §,—5 and 1 /^ = 1 
for a = i, — I, respectively, in which the chiral operator 
is defined as Cch = iCpCx = ® R. The mean-field 

Hamiltonian H{k\\,z) is in the Dili class satisfying the 
particle-hole and TR symmetries, and it transforms as 
CchH{k\\,z)C~^ = -H{k\\,z). Thus Cch is a symmetry 
only for zero modes. For a nonzero mode |'0n) and its 
chiral partner |' 0 fi) = Cch\'4’n)^ their energies are oppo¬ 
site to each other, i.e., 6 ^ = —e„. If a perturbation 5H 
remains in the Dili class, then Cch5HC~f^ = —6H. SH 
can only mix two zero modes with opposite chiral indices 
because {0aj\SH\0pj) = (—and 
it is nonzero only if Va ^ vp. 

As moving away from fey = 0, the zero modes evolve 
to the midgap states developing energy dispersions. At 
fcy fc/, these midgap states can be solved by using the 
k ■ p perturbation theory within the subspace spanned 
by the zero modes |0aj) at fcy = 0. By setting SH = 
H (A:y ,z) — H (0, z), the effective Hamiltonian to the linear 
order of fcy is 

0 —ik- 0 0{k^) \ 

ik+ 0 —2ik- 0 | ,.. . 

0 2ik+ 0 -ik. ’ 

0(fc+) 0 ik^ 0 / 

where k± = kx A iky. The matrix elements in the same 
chiral sector are exactly zero, and the elements at the 
order of 0{k\) are neglected. The solutions consist 
two sets of 2D surface Dirac cone spectra represented 
by E'^^\k\\) = ±Va{b)k\\. The velocities are solved as 
^ ^)- develop a systematic 

method beyond the k ■ p theory to solve the midgap spec¬ 
tra for all the range of fcy as presented in Appendix (HI 
and the results are plotted in Fig. [5] (a). In addition 
to the Dirac cones, there also exists an additional zero 
energy ring not captured by Eq. [TI] which is located at 
k/kf = as analyzed in Appendix [E] 

Now we move to other unconventional isotropic pair¬ 
ings of spin -1 fermions in the p and d-wave channels. 
The p-wave triplet one is topologically non-trivial, and 
the analysis can be performed in the same way as above. 
The pairing matrix is Ap(fc) = = 

Ap-^Kp(k)R, where is the rank -1 spin tensor, dp = 
Ap(^)Yii/(fc), and Kp(k) = k ■ S is just the helicity op¬ 
erator. The quasi-particle spectra are fully gapped as 
E\(k) = ^e^{k) + \ Ap\^{-^YX^, and the topological in¬ 
dex of this pairing can be evaluated based on Eq. [9] by 
replacing the eigenvalues of Kf(k) with those of Kp(k). 
The contributions from two helicity pairs of A = ±| and 
±5 add up leading to a high value = 4. In compar¬ 
ison, the topological index of the ^He-B phase is only 1, 


HLdih) = 


Akf 
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and thus their topological sectors are different in spite of 
the same pairing symmetry. 

The surface spectra of the isotropic p-wave pairing 
with spin-1 fermions are interesting: They exhibit a cu¬ 
bic Dirac cone in addition to a linear one. Consider the 
same planar boundary configuration as before, similarly 
for each spin component a there exists one zero mode at 
fc|| = 0 labeled by |0 q,_p). Again we perform the k-p anal¬ 
ysis at fc|| -C /c/ in the subspace spanned by | 0 a_p) with 
respect to SH = iJ(fc||,z) — H{0,z). The chiral eigen¬ 
value of | 0 q,p) is (—)'^“ = sgn(a), which leads to a dif¬ 
ferent structure of effective Hamiltonian from that of the 
/-wave one. Only |0±i p) can be directly coupled by 6H, 
which leads to a linear Dirac cone. In contrast, the pair 
of states | 0 j _3 p) are not directly coupled, rather |03_p) 
and | 0 _i p) are coupled through the 2 nd order perturba¬ 
tion theory, and so do | 0_3 p) and |0i p). Consequently, 
| 0 j _3 p) are coupled at the order of {SH)^ developing a 
cubic Dirac cone as shown in Appendix [F] The above 
analysis is confirmed by the solution based on the non- 
perturbative method in Appendix |h 1 as plotted in Fig. 
[ 21 (b). 

In contrast, the d-wave spin quintet isotropic pairing 
of spin-1 fermions are topologically trivial. By imitat¬ 
ing the analyses above, we replace Kf{k) with Kd(k) = 
2(k ■ Sy — I/ 4 . Different from the kernels Kp and Kf 
in odd partial-wave channels, K^s eigenvalues are even 
with respect to the helicity index, i.e., such 

that vanishes. This result agrees with that 3D TR 
invariant topological superconductors should be parity 
odd as shown in Ref. The explicit calculation of 

the surface spectra in Appendix |G] confirms this point 
showing the absence of zero modes. 

The above analysis can be straightforwardly applied 
to multi-component fermion systems with a general spin 
value s = The spin-tensors at the order of I are de¬ 

noted as 5*™ with 0 < I <2S and —l<m<l. For each 
partial-wave channel 0 < Z < 2S', there exists an isotropic 
pairing with the pairing matrix A(fc) = Ai{-^yKi{k)R 

in which Ki = U{k)S’'^U^ (k), whose topological index 
fVu,(Z) is determined by the sign pattern of the elements 
of the diagonal matrix For even and odd values of 
Z, respectively, and thus vanishes 

when Z is even, while for odd values of Z, 

ZV^(Z) = ^ 2 Asgn(^i'i), ( 12 ) 

A>0 

in which = (—)“+5 (Sa, S—a\SS] 10) up to an overall 
factor. The largest value of is reached for the p-wave 
case: Since oc Sz, contributions from all the com¬ 
ponents add together leading to iVy, = N'^. The ^He-B 
phase of spin-i fermions and the isotropic p-wave pairing 
with spin-1 fermions are two examples. As for the sur¬ 
face zero modes jOa^p) at fc|| = 0, their chiral indices equal 
sgn(a). As a result, similar to the spin-| case, when 


performing the k ■ p analysis for midgap states within 
the subspace spanned by |0a_p), only |0j_i p) are directly 
coupled leading to a linear Dirac cone, and other pairs 
of |0±Q,p) are indirectly coupled at the order of 
leading to high order Dirac cones. 

Multi-component fermion systems are not rare in na¬ 
ture. In solid state systems, many materials are orbital- 
active including semiconductors, transition metal oxides, 
and heavy fermion systems. Due to spin-orbit coupling, 
their band structures are denoted by electron total an¬ 
gular momentum j and in many situations J > 5 . For 
example, in the hole-doped semiconductors, the valence 
band carries j = | as described by the Luttinger model 
( 4 ^ . Superconductivity has been discovered in these sys¬ 
tems including hole-doped diamond and Germanium |46l - 
[ 4 ^ . Although in these materials, the Cooper pairings are 
mostly of the conventional s-wave symmetry arising from 
the electron-phonon interaction, it is natural to further 
consider unconventional pairing states in systems with 
similar band structures but stronger correlation effects. 
The p-wave pairing based on the Luttinger model has 
been studied in Ref. [4l|. In ultra-cold atom systems, 
many alkali and alkaline-earth fermions often carry large 
hyperfine spin values F > ^, and thus their Cooper 
pair spin structures are enriched taking values from 0 
to 2F not just singlet and triplet as in the spin-i case 

[Ulsllil. 

In multi-component solid state systems, there often 
exists spin-orbit coupling. For example, the Luttinger 
model describing hole-doped semi-conductors [i^, con¬ 
tains an isotropic spin-orbit coupling FI^o = 72 Zc^(fc • S)'^. 
Since iLgo is diagonalized in the helicity eigenbasis, we 
only need to update the kinetic energy with ek\ = 
Cfc + in the mean-field analysis, which satisfies 

CfcA = Cfc.-A, and the pairing structure described by Eq. 
[ 6 ] is not affected. The topological properties are the same 
as analyzed before because the index formula Eq. [9] re¬ 
mains valid and the surface mid-gap state calculation 
can be performed qualitatively similarly. Nevertheless, 
the symmetry breaking pattern is changed. The relative 
spin-orbit symmetry is already explicitly broken by the 
Flso- The spin-orbit coupled Goldstone modes in ^He- 
B become gapped pseudo-Goldstone modes with the gap 
proportional to the spin-orbit coupling strength 72 . 

In summary, we have found that multi-component 
fermion systems can support a class of exotic isotropic 
pairing states analogous to the ^He-B phase with uncon¬ 
ventional pairing symmetries and non-trivial topological 
structures. High-rank spin tensors are entangled with 
orbital partial-waves at the same order to form isotropic 
gap functions. Eor the spin-| case, the odd partial-wave 
channel pairings carry topological indices 2 and 4 for 
the / and p-wave pairings, respectively, while the cZ-wave 
channel pairing is topologically trivial. The surface Dirac 
cones of mid-gap modes are solved analytically which ex¬ 
hibit two linear Dirac cones in the /-wave case, and the 
coexistence of linear and cubic Dirac cones in the p-wave 
case. Generalizations to systems with even more fermion 
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components can be performed straightforwardly. This 
work provides an important guidance to search for novel 
non-trivial topological pairing states in both condensed 
matter and ultra-cold atom systems. 
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Note added. After the submission of this manuscript, 
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has been reported in the rare earth-based half-Heusler 
superconductors 


Appendix A: Spherical harmonic functions and 
high-rank Spin tensor operators 

In this section, we present spherical harmonic functions 
in momentum space and high-rank spin tensors. 

For convenience in the main text, we normalize the 
spherical harmonic functions Yim{k) defined on the Fermi 
surface satisfying 


^ = 1 . 


(Al) 


m— — l 


This normalization differs from the usual o ne o f 
/ = 1 only by an overall factor 

More explicitly, for the p-wave case, they are defined as 

kYi±i{k) = =F-^fc±, kYio{k) = kz, (A2) 


For the d-wave case, they are defined as 

/3, 


k Y2±2{k) — \ k Y2±i{k) — —k±kz, 


k^Y2oCk) = ^(3fc2-fc2). 

For the /-wave case, they are defined as 


(A3) 


k^Y3±sCk) = T^kl, k^Y3±2ik) = ^k^kz, 
k^Y3±iik) = ±^fc±(fc"-5fc^), 
k^Ysoik) = -^{3e-5kl)kz. 


(A4) 


All of them are homogeneous polynomials of momentum 
components kx, ky and kz- 

The spin-1 matrices are defined in the standard way 


as 




/O 

Vs 0 

0 \ 

5+ 


0 

0 

2 

0 

— 

0 

0 

0 

V3 



VO 

0 

0 

0 / 

5_ 

— 






/ ^ 
2 

0 

0 

0 \ 

5, 


0 

1 

0 

0 

— 

0 

2 

0 

1 

2 

0 



^0 

0 

0 

-1) 


(A5) 


in which S± = Sx^- iSy. The general rank-A: spin tensors 
Sjm satisfy 

[S-,Sjm] = V(j + m){j - m + l)Sjm-i- (A6) 

Based on these relations, we can build up spin tensors. 
For example, the rank-1 tensors are defined as 

Sn = -^S+, Sio = Sz, 5i_i = ^5_. (A7) 

The rank-2 tensors are spin quadrupole operators de¬ 
fined as 


^22 


S 2 I — 2[‘5'-)'5'22] = 


/O 0 1 0 
0 0 0 1 
0 0 0 0 
Vo 0 0 0 

/O -1 0 0 
0 0 0 0 
0 0 0 1 
\ 0 0 0 0 . 


520 ^ 


/ 1 0 0 0 
0-100 
0 0-10 
Vo 0 0 1 


52-1 — "^^[5-, 52o] — —521 


1 , 


52-2 = ^[5_,52-l] = 5|2- 


(A8) 


This set of tensors can be organized into the Dirac F 
matrices through the relations of 

Fi = —*(522 — 52-2)5 r5 = 522 + 52-2 
+2 = —521 + 52-1, +3 = *(521 + 52-i) 

+4 = V2520, (A9) 


which satisfy the anti-commutation relation 

papb pfcpa ^ 2 Sab. 


(AlO) 


The rank-3 spin tensors also called spin-octupole 
operators, are constructed as follows 
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Q c3 

~ 


^3,-1 

83,-3 


/O 0 0 1 
0 0 0 0 
0 0 0 0 
Vo 0 0 0 


^31 = 




/O 1 0 OV 

0 0 -^3 0 
0 0 0 1 
\ 0 0 0 0 J 


1 


2 V 3 
1 

To 


[S.,S 3 o] = -Sl, 


83,-2] = -S} 


33 




830 = 


2 V 3 


[5-, 531] = 


2 V 5 


/O 0 1 0 
0 0 0 -1 
0 0 0 0 
Vo 0 0 0 

/ 1 0 0 
0-3 0 
0 0 3 


V 0 0 



83.-2 = 


Vio 


[8., 83,-1] = 8: 


32 7 


r 


(All) 


Appendix B: Pairing matrices for the isotropic p, d, 
and /-wave 


matrix kernel Yi^{k)8''™ can be explicitly represented in 
isotropic forms as 


By using the spherical harmonic functions Yim(k) and 
spin-tensors, we can construct the pairing matrices for 
the isotropic pairings for the spin-| fermions in the p, d, 
and /-wave channels, respectively. The pairing matrix 
^a/sik) in momentum space can be represented as 


k- 8 , (1 = 1 ) 

YL{k)8^^ ={ • Sr - (/ = 2) 

fik-Sr-^^k\k-S), (/ = 3) 


(B2) 


Kpik) = qA i-rY,mik)8^^R 

= ciA Y;^{k)8^^R, (Bl) 

in which I = 1,2,3 represent p, d, and /-wave pairings, 
respectively, while c; is an overall constant factor. These 
pairing structures are isotropic in analogy to the ^He-B 
phase: the pairing orbital angular momenta are I, and 
the pairing spins are also I, such that they add together 
into the channel of total angular momentum J = 0. The 


More explicitly, the pairing matrix in momentum space 
can be expressed as follows. In the p-wave case (ci = 1), 
it is 


^k^ 

2 ^+ 
0 



/ 

0 

0 

1 

1 

A 


0 

k_ 

1 '-i|C^ 
1 

kf 


V3u 

2 

--k 

2 '^z 

-k+ 


1 

Ik 

2^z 

2 

0 


0 J 


(B3) 


J 


In the d-wave case (c 2 = 2-\/2), it reads 


^Uk) = 


kj 


0 

— 

2i/3k-kz 
\ k^ — 3k^ 


•\/ 3 fc^ 

—2y/8k-kz 

3kl — k"^ 

\ 

0 

1 

CO 

2^/3k+kz 


— 8k1 

0 

-J3k\ 


—2\/8k+kz 


0 

/ 


and in the /-wave case (03 = — it becomes 


^ 3V3A 

8k} 


( 


W- 


—un,_n,z —\/8k-{k8 — ^k1) V3{3k‘^ — 5k^)kz k+{k'^ — 5k'^) 
—— 5fc^) '/8{‘ik8 — bk1)kz \/3k+(k'^ — 5k^) 

V k+iP-bk}) 


—k-{k‘^ — bk‘^) -^{3k‘^ — 5k‘^)kz \ 

/ 


—^k'}_kz 


—bk\kz 


_At3 


The general mean-field Hamiltonian in coordinate space is represented as 


H = 


d^f'^^{r} 


Cl-, 




«'(r), 


(B4) 


(B5) 


(B6) 
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in which '1''^ (r) = (cjj(r), C/ 3 (r)) is the Nambu spinor; —zV 
replaces k in the expressions of k!‘Yim{k). 


Appendix C: The symmetry properties 


Here we use the isotropic /-wave pairing state as an ex¬ 
ample to illustrate the symmetry properties of this class 
of pairings. The Hamiltonian Eq. 1 (main text) pro¬ 
cesses both rotation symmetries in the orbital and spin 
channels, be., S'Ol(3) x S'Os( 3), while in the isotropic 
pairing state characterized by Eq. 5 (main text) only the 
total angular momentum is conserved, z.e., the residue 
symmetry is SOj{2>). The general configuration of the 
d-tensor can be expressed as 

d^(fc) = d‘'{R-^k) = Dl=f{R)d'''{k), (Cl) 

where R is an arbitrary SO(3) rotation and DIj^,{R) is the 
rotation H-matrix. The relative SO symmetry is spon¬ 
taneously broken similar to the case of ^He-B, and here 
it is realized in a high representation of angular momen¬ 
tum. Combining with the 17(1) gauge symmetry breaking 
in the paired superfiuid state, the Goldstone manifold is 
[SOl{5) O SOs{5) O 17(1)]/S'Oj(3) = S'0(3) O U{1). Ac¬ 
cordingly there exist four branches of Goldstone modes, 
including one branch of phonon mode and three branches 
of relative spin-orbit modes. 


spending eigenvalue A is defined as 


q\ = 


k^dilk 

Air 

k'^dQk 


FeMk) 


j (VeC/V^t/t - V^UVeU^) 


\X 


= A, 


(D3) 


and w\ is the winding number of the angular 9\{k) along 
the radial direction of k as 


wx 


Jo 


(*VAAt)^^ 


sgn(^A) (m > 0 ), 

0 (Ai< 0 / 

(D4) 


Consequently, we arrive at 


/EA^sgn(/A) (m>0). 

\0 (/i< 0 ). 


(D5) 


Appendix E: The surface modes of the /-wave 
isotropic pairing 

In this part, we study the gapless surface states of the 
/-wave septet pairing. 

We study a boundary imposed at z = 0 with a spatial 

hk^ 

dependent chemical potential d{z): > 0 at z < 

0 and /Tiz < 0 at z > 0. Eor simplicity, we consider the 
case of /iL and finally take the limit of —>■ oo, 

i.e., at z > 0 is the vacuum. 


Appendix D: Calculation of the bulk topological 
index 


1. The /-wave zero modes at fc|| = 0 


In this section, we calculate the topological index of 
various pairing states. According to the definition of 
Q(k) = W{k)A{k)U{k) in which U and only depend 
on the direction of k, while A{k) only depends on the 
magnitude of fc, we have 


To warm up, we first consider the case of k\\ = 0 in 
which Sz remains a good quantum number, and the zero 
modes described by different Sz eigenvalues decouple. 
The B-deG equation of the zero mode with S'^-eigenvalue 
a becomes 


VfcQ = U^VkA{k)U, 

VeQ = VeC/U{7-f C/UVeC/, 

= V^C/UC/+ C/UV^C/, (Dl) 

in which Vfc = fc • V, Vg = • V, ■ V. Substi¬ 

tuting the above equations into Eq. 8 in the main text, 
after simplification, we arrive at 

J d^kTiiVeUV^U^ -V^UVgU^)VkAA\ 

= '^q\wx, (D2) 

A 

in which qx is the monopole charge associated to the 
Berry curvature of the helicity eigenstate; the corre- 


0 , 

(El) 

in which Aj_3 = |A and Aj_i = |A. The boundary 
condition is that 

uo{z) —>■ 0, vo{z) —>■ 0, (E2) 


2m dz^ 

;A„ h? d? ' 


d-' 

dz^ 


-I- 


k% dz^ 


2m dz^ 


■Ai(z) 


<(z) 

V-a{z) 


as z —>■ ±oo. 

Eq. 10 (main text) is invariant under the operation 


of 


uo 

Vo 


—>■ iT2 


Uo 

Vo 


in which T 2 acts in the Nambu 


space, thus we can set vq = itiuo- As it will be clear later 
that the solution actually satisfies vo(z) = —iuo(z), the 
other one with vo(z) = iuo(z) corresponds to the case 



that the system lies at z > 0 and the vacuum is at z < 0 . 
Then the equation becomes 

f fi^ cP , , Aq, dp \ 

If 5? j = "^ <“> 

We try the solution uq{z) ^ at z < 0, and, at 
z > 0 , and thus Re/Sn > 0 and Ro/Sr < 0 , respectively. 

At z < 0, Pl satisfies the cubic equation with real 
coefficients 


kf 


3 




(E4) 


which has a pair of conjugate complex roots and one 
real root. We only consider the weak pairing limit that 
— <C 1. The solutions correct to the linear order of 

^ < 1 are 


Pl 

kf 


±i + 


1,2 


2 e/’ 




(E5) 


and thus only (§^)i ,2 can be kept. Similarly, at z > 0, 
in the case of |/j,r| ^ e/, there exists a pair of complex 
conjugate roots and one real root for /3 r as 


pR 

kf 


1 , 


(E 6 ) 


in which c = (l/iRl/Aa) ^. 

Because Eq. 10 (main text) is a 3rd order differential 
equation, all of uo(z), -^uo{z), ^uo(z) need to be con¬ 
tinuous at the boundary z = 0. For this purpose, we 
construct the following solution 

uo{z) = l (z < 0 ) ^ 

[ Ar sin(-^cfc/z-I-(/)R)e 2 ^/^^ (z > 0) 

in which the four parameters Al{r) 4’l(r) sic suffi¬ 
cient to match three continuous conditions. In the case 
of c —>■ - 1 - 00 , the results can be simplified as 


Ar 

Al 


0 , (pR = 

1 

csin(| - 4>r) 


(E 8 ) 


which shows that we can simply set ul{z) vanishing at 
z = 0 . 

To summarize, we have solved 

w/.a( 2 :) =-;^^e^“^sinfc/z (E9) 


with / 3±3 = and = 3A±, = 

3 |A| kf ■ 


2. The k ■ p perturbation theory for midgap states 


The effective Hamiltonian for the midgap states on the 
surface of the /-wave isotropic pairing is presented in Eq. 
11 (main text). The spectra consist of two gapless Dirac 
cones denoted as a and b, respectively, as shown in the 
main text. The corresponding eigenfunctions are solved 
as 


r±{h) 


V'±(fc||) 



( ipjg \ 


1 

— i^Js. 


i 

—xe 2 



±ixe'^~ 

5 





/ ±ixe \ 

1 

• 1 

g '^ 2 Vk 





/ 


(ElO) 


in which x = y/2+1; N = 2\/2 + \/2; (j)k is the azimuthal 
angle of fey. 

The eigen-solutions are parity eigenstates of 

the little group symmetry of the reflection cr„(fc||), which 
is defined with respect to the vertical plane passing k\\ 
and the z-axis z. The operation cr«(fc||) can be decoupled 
as a combined operation of inversion and rotation as 


tT„(fc||) = iIRp{'K) 


/ 0 0 0 \ 

0 0 0 

0 - 6 *'^'“ 0 0 

\ 0 0 0 / 

(Ell) 


in which / is the inversion operation; (j)k is the azimuthal 
angle of fey; k' is an in-plane momentum perpendicular 

to k\\ and (tt) is rotation around k\\ at the angle of tt; 
the factor of i is to make (t„ an Hermitian operator with 
eigenvalues ±1. It is easy to check that 

cr^(fcy)V'±(fc||) = ±V'±(fc||), 
cr^(fcy)V'±(fc||) = TV'±(fc||), (E 12 ) 


respectively. 


3. The surface zero energy ring states 

Here we present the eigenfunctions of the zero en¬ 
ergy ring of the midgap surface states of the isotropic 
/-wave pairing state, which is located at fc!? = ^kf. The 
method of solution can be referred to SM[H1 The zero en¬ 
ergy states at fcy = fcjj(cos^fc,sinc/ife) are two-fold degen¬ 
erate, whose creation operators are denoted as 71 , 2 (^ 11 ), 
respectively. They are explicitly expressed below as 
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7l(fc||)= r dz ^fe*(^+^)4(fc||,^) + (-)“+^e-*(^+T)c„(-A||,^) 


' —oo 
rO 




72 (^ 11 )= / dz '^\{-)°‘ 2 e *(2 + I)ct^(fc||,z)+ e *(2 + I)c„(-fc||,z) e *“'^'=M_a(z) 

00 


(E13) 


in which (j) is the phase of the pairing amplitude A; 
a = zb|,zb^ as the eigenvalue of 5^; the envelope wave- 
functions Ua{z) are 


Uslz) = 

2 ^ ' 

Ui(z) = 

U_i(z) = 
2 ^ ' 

U_3(z) = 


1 - 


3cos(^)(e^^^ 


Ns 

1 


sin(^)(3e'5i" + 5e'5="), 


= -\/3cos(%^)(e^^^ — 
N 1 2 


N_ 




(E14) 

in which (32 = 3^i = and NaS are the overall 

normalization factors whose expressions are complicated 
and will not be presented. 

Since 7^ 2 (^ 11 ) represent the zero energy modes, again 
they are chiral eigen-modes satisfying 


Cchli{k\\)Cj = 7i(fc||), 

a/.72(fc||)C-^' = -72(fc||). 


(E15) 


Nevertheless, they are not parity eigen-modes, and they 
transform into each other under the parity operation de¬ 
fined in Eq. lElll as 


CT«(fc||)7i(*||)o-« ^(fc||) = 72(*||)- 


(E16) 


The equation parallel to Eq. IE3I is 

AA _ f i A 

2m dz"^ ^ kf dz 


uoiz) = 0. (El) 


in which Aj _3 = |A, Aj_i = —^A, and vo{z) = iuo{z) 
for a = ±|, vo{z) = —iuo(z) for a = Again we set 
uo(z) ^ at z < 0 , and, at z > 0 , with Rc/Sl > 0 
and Re/^R <0. At z < 0, /3 l satisfies the equation 


kf 


e/ 




which in the limit — << 1 has the solutions 




f 7 1,2 


2 e/ 


(F2) 


(F3) 


At z > 0, in the limit of \fLR\ » e/, /3 _r has two solutions 
as 


Pr,i,2 


± 


V £/ 


(F4) 


and the negative one is kept to match boundary condition 
at z —>■ -boo. 

Eq. m is a second order differential equation, and 
thus Uo{z) and -^Uq{z) need to be continuous at z = 0 . 
Similar to the /-wave case, we arrive at 


Up,a{z) = sin/c/z, 

with and jA. 


(F5) 


|A„| kf- 


Appendix F: The surface states of the p-wave 
isotropic pairing 

In this section, we consider the surface states of the p- 
wave case under the same planar boundary configuration 
as that in the /-wave case. 


1. The p-wave zero modes at fc|| = 0 

There also exists one zero mode at fey = 0 for each spin 
component in the p-wave case, whose spatial wavefunc- 
tions will be solved below. 


2. The k ■ p perturbation theory 


The chiral indices for the zero modes |0 q,)p at ku = 0 


are 1,1,—1,-1 for a = respectively. We 

can use these zero modes as the bases to construct the 
effective Hamiltonian for low energy midgap states at 
fc|| <C kf with respect to 6H = i7(fc||,z) — H{0,z). As 
constrained by the surface symmetry C„oo and the chiral 
symmetry of the zero modes, the effective Hamiltonian is 


/ 


HLM = 


0 
0 

ck\ 

\0{kl\ 


0 ckl 0{kl_) 

0 —ifc_ 

ik+ 0 0 

ck\ 0 0 


(F6) 
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in which the terms proportional to arise from the 3rd 
order perturbation theory and are neglected. The terms 
proportional to are due to the 2 nd order perturba¬ 
tion theory involving gapped bulk states as intermediate 
states. Under a suitable phase convention c is a real co¬ 
efficient at the order of l//c/, whose concrete value is not 
important. At the leading order, the two components 
with a = ±i form a linear Dirac cone, while the other 
two with a = ±1 are dispersionless. Nevertheless, the 
latter are coupled indirectly through coupling with the 
former and develop a cubic Dirac cone. 


Appendix G: The isotropic d-wave pairing 


The isotropic d-wave pairing with spin-| fermions is 
actually topologically trivial. In this part, we explicitly 
check this point from the boundary spectra. The bound¬ 
ary configuration is the same as before, and we will show 
the absence of the zero modes at fey = 0 . 

Similar to the /- and p-wave cases, the equation deter¬ 
mining zero modes is invariant under zt 2 operation. Let 
vo{z) = iriuo{z) (rj = ± 1 ), we obtain 

(-^ - uo{z) - = 0, (Gl) 

in which As = Ai = 2A, A_3 = A_i = —2A. Ex- 

2 2 2 2 

pressing uo{z) ^ at z < 0 , and, at z > 0 , Pl 

and Pr are solved as 



s ± 

kf 



s ± 

kf 



A„ 
1 — 
e/ 


£/ V 2*^ €f 


(G2) 


where |Aa| <C e/ <C IphI is assumed. 

Since Eq. EH is a 2 nd order differential equation, both 
uo{z) and -^Uq{z) need to be continuous at z = 0. Re¬ 
gardless of value of ry, there is only one Pl with a pos¬ 
itive real part, and one Pn with a negative real part. 
The boundary conditions at z = 0 are two linear ho¬ 
mogeneous equations of the undetermined coefficients of 
wavefunctions. Generally speaking, there is only zero so¬ 
lution, which demonstrates the absence of zero modes for 
d-wave isotropic pairing. 


Appendix H: The method for solving the midgap 
surface states 

In this section, we present a general method for solving 
surface states in the weak pairing limit away from the 
A:|| = 0 point. The isotropic /-wave pairing in the spin-| 
system is used as an example, and actually the method 
can be directly applied to other partial-wave channels 
and higher spins. 


1. Match boundary conditions 


Gonsider the same boundary configuration as stated 
before in Supp. Mat. m We denote $( 7 ^ the eigen- 
wavefunction with the eigen-energy E and the in-plane 
wavevector /cy = (kx^ky). The following trial solution 
will be used 


<I>(r) 


z^ikxX+ikyy z: < 0 

^R^P^z^ik^x+ikyy ^ > 0 


(HI) 


in which <I>^, are 8 -component column vectors. De¬ 
note Hl and Hji the Hamiltonians for z < 0 and z > 0 
with the corresponding chemical potentials and p-R, 
respectively. Substituting the trial wavefunction into the 
eigen-equations at z < 0 and z > 0 , respectively, the con¬ 
ditions for the existence of nonzero solutions of and 
are obtained as 

J det[HL{kx, ky, -iP^) - A] = 0, . . 

\ de%[HR{kx,ky, -iP^) -E]=0. 

Both solutions —ip^ and —iP^ appear in terms of com¬ 
plex conjugate pairs, since the determinant Eq. IH2I are 
real equations. Consequently, among the 24 solutions of 
P^, there are 12 solutions with positive real parts and 12 
with negative real parts, and so do P^’s. The midgap 
state needs to vanish at z —>■ ±oo, hence, it is in the form 
of 


$(r) = 

in which Re/3^" > 0 and Re/3?^ < 0 , 1 < j < 12 . 

The boundary conditions require the wavefunctions 
Eq. IH31 and their first and second order derivatives to be 
continuous at z = 0. We have a set of linear homogeneous 
equations that the coefficients and should obey. 
The conditions for the existence of nonzero solutions are 

[El Er\ 

det I Er Fr 1=0 (H4) 

\Gl Gr J 

in which Er^r, Fr^r and Gr^r are 8 x 12 rectangular 
matrices, and thus the total dimension is 24 x 24. The 
above block rectangular matrices are expressed as 

Er{;j) = <!>^, EK(-,j)=$f, 

FrP.j) = pf<^f, Ffl(.,j) = /3f4>f, 

Gl{;j) = GR{;j) = {p^)Hf, (H5) 

in which (•,/) denotes j-th column of the corresponding 
matrix. Surface energies can be solved from this equa¬ 
tion. 

Actually the complicated determinant equation Eq. 
IH4I can be greatly simplified in the weak pairing limit 
as will be shown in Sect. IH 21 


z < 0 , 
z > 0 , 





11 


2. Equations of the midgap state energy 

In the half space at z < 0, we rewrite the eigen-solution 
<i)(r^ in Eq. IHSI as 

12 

i=i 

in which kzj = (1 < j < 12), hence, Imk^j < 0 

such that $(r) vanishes at z —>■ —oo. 

It can be shown that the twelve can be classified 

into two groups. In one group, their real parts are very 
close to ~ as 

kf^^ = ±^k}-kl-^S,t^, (H7) 

in which m = ±|, ±^, and are corresponding eigen¬ 
vectors. The remaining four k^'s represent fast decaying 
modes in the weak pairing limit, which are proportional 
to —i^kf. It can also be proved that at the leading order 
of the wavefunction at z < 0 is represented as 

$(r) = ^ (H8) 

m 

I 


in which the fast decaying mode contributions are ne¬ 
glected. 

It can be shown that in the limits of A <C e/ -C |/r_R|, 
the boundary conditions can be further simplified as 
<I>(r^ = 0 at the boundary of z = 0. The detailed proof is 
rather complicated but straightforward, and will be pre¬ 
sented elsewhere. This great simplification reduces the 
equation determining surface energies from the original 
24 X 24 determinant condition of Eq. IH4l to the following 
one of 8 X 8, 

det(^'+,4'-) =0, (H9) 

in which '1'=*= are 8x4 matrices defined by 4/^ (•, m-f |) = 
4)^^ (m = ±|, ±1). Furthermore, due to the SU{2) bulk 

rotation symmetry and the reflection symmetry ay{k\\) 
defined in Eq. lEIIl Eq. IH9I can be further simplified to 
two 4x4 determinant equations. 

The surface midgap energies are smaller than A, we 
express E = eA. To solve kf„i kf = 

— fey — i^^A is plugged into the eigen-equation. 
Keeping only the leading order of A, we obtain 


{Kik,,ky,±^kj-kpR)^ ±2iA^k} - k^C^h ) 


(HIO) 


in which the subscript m is dropped for simplicity. Since 
K{kx,ky,—idz) already contains a prefactor of A, the 

—idz can be substituted by iy/fcj — fcy without induc¬ 
ing higher order error. Denote U±{k) as the rotation 
matrices associated with the operations rotating ±z to 
the direction of k, the pairing part can be represented as 

K{±kfz)R = U^Hk±) K{k±)R C/±(fc±), (Hll) 


in which k± = {kx, ky,±^k‘j — fcy). Applying such rota¬ 
tion operations to the eigen-equation, we obtain 


/ T2zA^fc^ - k^i^h K(±kfz)R 
(K(±A;/z)i?)t ±2iA^k) - 


= £A$^^ 


(H12) 


in which 


and W is defined as 


= IT±^(fc±)4> 


hL± 


(H13) 


VE±(fc±)= 


U±(k^) 


Ul(k^ 


(H14) 
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Hence it is sufficient to solve to arrive at , which 
satisfy a simple equation where the pairing matrix is in 
±z direction. 

For notational convenience we define the column vec¬ 
tors representing particle and hole states as 

P| = (1000,0000)^, pi = (01000,000)^, 

p_i = (0010,0000)"^,p_3 = (0001,0000)^, 


and 

hs = ( 0000 , 1000 )'^, = ( 0100 , 0100 )^, 

h_i = ( 0010 , 0010 )'^, /i_| = ( 0001 , 0001 )'^. 

The solutions to and are summarized as follows 
(vectors un-normalized), 



Correspondingly, the determinant equation for the 
eigen-energies becomes 

det [lF+(fc+)l'^+,lF_(fc_)l'^- =0, (H16) 

in which are 8x4 matrices defined by ^^(-,771 -|- 
|) = (^ = As mentioned before, in this 

way, the original 24 x 24 matrix determinant equation is 
reduced to an 8 x 8 one. 

Further using the reflection symmetry, the above 8 x 8 
matrix can be further decomposed into two 4x4 ones. 
Without loss of generality, we only consider fey along the 
x-axis, i.e., /cy = (fcy,0), and results for other values of 
k\\ can be obtained by applying rotations around the z- 
axis. The reflection operator with respect to the vertical 
xz-plane which we denote by ayx, is given in the particle- 
hole 8 -dimensional space as 

= ( 1 ” « ) ^ 

The vectors can be recombined into even and odd 
eigenvectors of ayx defined as 

’2 2 2 

2 2 2 

^’2 2 2 

(H18) 

^’2 2 2 


in which the subscripts “e” and “o” denote even and odd 
parity eigenvalues 1 and —1 of a^x, respectively. 

For fcy = (fcy,0), U±{k±) are rotations around the y- 
axis, which commutes with ayx- Applying a basis trans¬ 
formation P which separates the even and odd parity 
eigen-spaces of ayx, we have 

= (1”)' ^"b*, = (.jy.(Hi9) 

where y = §, 5 , and then 

(H20) 

In this set of basis, we obtain the following two 4x4 
determinant equations for the even and odd sectors of 
Uyx, respectively, as 

det (l7+,e^^,C/_,e^e") = «, 

det ([/+,„§+, =0, (H21) 


in which U±^e, 



U±^o are 4x4 matrices, and = 
= (<('q 3 , 1 ) are 4x2 ones. The sur¬ 


face midgap state spectra displayed in the main text are 


solved from this set of equations which are fourth order 
algebraic equations of e^. 
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